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Abstract. The non-gaussianity of processes observed in financial markets and rela- 
tively good performance of gaussian models can be reconciled by replacing the Brownian 
motion with Levy processes whose Levy densities decay as exp(— X\x\) or faster, where 
A > is large. This leads to asymptotic pricing models. The leading term, Pq, is the 
price in the Gaussian model with the same instantaneous drift and variance. The first 
correction term depends on the instantaneous moments of order up to three, that is, 
the skewness is taken into account, the next term depends on moments of order four 
(kurtosis) as well, etc. In empirical studies, the asymptotic formula can be applied 
without explicit specification of the underlying process: it suffices to assume that the 
instantaneous moments of order greater than two are small w.r.t. moments of order 
one and two, and use empirical data on moments of order up to three or four. As an 
application, the bond pricing problem in the non-Gaussian quadratic term structure 
model is solved. 

For pricing of options near expiry, a different set of asymptotic formulas is developed; 
they require more detailed specification of the process, especially of its jump part. The 
leading terms of these formulas depends on the jump part of the process only, so that 
they can be used in empirical studies to identify the jump characteristics of the process. 
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I. Introduction 

To account for fat tails, skewness and excessive kurtosis of empirical probability dis- 
tributions of returns in real Financial Markets, it has become increasingly popular to 
model the dynamics of market factors as a Levy process. Levy models are more realistic 
than Gaussian ones but the latter are much more tractable. Indeed, in the Gaussian 
framework, explicit pricing formulas are known for a wide range of options and other 
contingent claims both without and with early exercise features, whereas in the Levy 
models, most of the pricing formulas have been obtained for contingent claims of the 
European type, with the deterministic life-span. There are some explicit analytic re- 
sults for options with early exercise features: see Boyarchenko and Levendorskii (2000, 
2001, 2002a, b), Mordecki (2002) and the bibliography therein for pricing of perpetual 
American options, and Boyarchenko and Levendorskii (2002b, c) for pricing of barrier 
options and first touch digitals. However, the pricing formulas are complicated and diffi- 
cult for numerical implementation except for a rather special case of pricing of perpetual 
American options under exponential jump-diffusions or spectrally one-sided processes. 

Another obstacle for non-Gaussian modelling arises when one considers more general 
Markov processes. The explicit pricing formulas in affine term structure models and 
certain Levy-driven Ornstein-Uhlenbeck models are known in the case of contingent 
claims with the deterministic life span only - see Duffie et al. (2000, 2002), Chacko 
and Das (2002), and Barndorff-Nielsen and Shephard (2001b), Barndorff-Nielsen et al 
(2002), respectively; for non-Gaussian variants of the HJM- model, see Eberlein and 
Raible (1999). In the general case, the dependance on the state variable does not allow 
one to obtain explicit analytical answers. 

The following observation helps to obtain efficient approximate solutions. As Barndorff- 
Nielsen and Levendorskii (2001) notice, typically, a good fit to the data can be achieved 
with Levy processes whose Levy densities decay as exp(— \\x\) or faster, where A > (the 
steepness parameter of the exponential Levy process) is large. They used this property 
to derive an asymptotic pricing formula for European options under certain class of Feller 
processes. The same observation was used in Boyarchenko and Levendorskii (2002a,b,d) 
and Kudryavzev and Levendorskii (2002) to derive efficient approximate formulas for 
perpetual American and Bermudan options, and first-touch-digitals, respectively. 

It was shown in Boyarchenko and Levendorskii (2002a, b) that the simple approximate 
formula is of the same form as the corresponding formula in a Gaussian model even when 
the underlying Levy process has no Gaussian component. It can be shown that the 
leading term of the approximate pricing formula in Barndorff-Nielsen and Levendorskii 
(2001) can also be written as the pricing formula in a Gaussian model. These observations 
can serve as an analytical explanation of relatively good performance of Gaussian models 
in apparently non-Gaussian situations. Thus, as far as pricing formulas are concerned, 
Levy processes with large steepness parameters behave almost as the Brownian motion, 
and Feller processes with large steepness parameters considered in Barndorff-Nielsen and 
Levendorskii (2001) behave almost as Gaussian diffusions. It seems reasonable to use the 
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nomer pseudo-diffusions for Levy processes and more general Levy-like Feller processes 
with large steepness parameters . 

The modelling with pseudo-diffusions allows one to obtain an efficient approximation 
to the price; in some situations, the asymptotic expansion of the price can be obtained, 
of the form 

(1.1) P(x, t) ~ P (x, 0(1 + X^Piix, t) + X~ 2 P 2 (x, t) + ■■■), 

where the leading term, P , is the price in the Gaussian model with the same instan- 
taneous drift and variance. The first correction term takes into account the moments 
of order three as well (skewness), the second correction term accounts for moments of 
order four, etc. Notice that though the leading term looks as the pricing formula in the 
Gaussian model, the "drift" and "variance-covariance matrix" used in the formula for 
the leading term are not the same as the ones of the Gaussian component of the process 
unless it is purely Gaussian. Indeed, a Levy process may have no diffusion component 
at all. 

The aim of the paper is to apply the approximate pricing approach to quadratic term 
structure models (QTSM) when the stochastic factor follows a mean-reverting pseudo- 
diffusion process of the simplest form (it is unlikely that in the QTSM model, an explicit 
pricing formula can be obtained unless the process process is Gaussian), and derive a 
pricing formula of the form (jl.lj) . For the discussion about advantages of the Gaussian 
QTSM model, see Ahn et al (2002, 2003) and Chen and Poor (2002). Cheng and Scaillet 
(2002) consider an affine-quadratic model, and allow for jumps but only in the dynamics 
of affine variables of the model. Notice that the use of jumps in QTSM models adds 
additional flexibility in joint modelling under the historic and a risk-neutral measures, 
and one may hope that the performance of QTSM models can be improved by introducing 
jumps. 1 Another improvement (and quite sizable one) is expected in pricing of out-of- 
the-money options near expiry, where the main contribution to the price comes from 
the jump part of the process. Near expiry, however, a different approximate formulas 
are needed, which use more detailed information about the jump part of the processes 
than the skewness and kurtosis. These formulas are similar to approximate formulas for 
out-of-the-money options on stocks developed in Levendorskii (2003), and can be derived 
by the same reasoning. 

1.1. Plan of the paper. In Section |21 we list families of exponential Levy processes 
used in empirical and theoretical studies of financial markets. In Section |21 we formulate 
the pricing problem for an interest rate derivative of the European type, and by using 
the Feynman-Kac theorem, reduce the pricing problem to the boundary problem for an 
integro-differential equation. We also explain the scheme of the asymptotic pricing. In 
Section we recall the solution of the bond pricing problem in the one- factor Gaussian 
case, and indicate the properties of the solution which are crucial for our asymptotic 
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method. In Section we demonstrate our method in the simplest case of the one- 
factor Levy model for the bond price, and present numerical examples. In Section El 
we consider possible specifications of the market price of risk, the generalization for the 
multi-factor case, derive approximate formulas for interest rate derivatives near expiry, 
and suggest a procedure of parameter fitting based on the asymptotic expansions. In 
Section [7J we summarize our results, and compare the Levy QTSM with multi-factor 
Gaussian QTSM. In the appendix, technical results are proven. 

2. Levy processes in financial modelling 

As early as in 1963, Mandelbrot suggested to use stable Levy processes. The modelling 
with stable Levy processes is not quite realistic since the tails of Levy stable distributions 
are too fat (polynomially decaying), whereas the tails of distributions of returns observed 
in real financial markets exhibit exponential decay. Moreover, the second moment of a 
Levy stable distribution is infinite (unless it is a Gaussian one). This contradicts the 
observed convergence to the Gaussian distribution over a longer time scale, and even 
worse, the underlying stock itself should have the infinite price under the stable Levy 
process, which makes the model inconsistent for pricing purposes. Starting with the 
beginning of the 90th, several families of Levy processes with probability distributions 
having exponentially decaying tails have been used to describe the behavior of stock 
prices in real financial markets: 

• Variance Gamma Processes (VGP) constructed and used by Madan and co- 
authors in a series of papers during 90th (see Madan et al. (1998) and the 
bibliography therein); 

• Hyperbolic Processes (HP) were constructed and used by Eberlein and co-authors 
(see Eberlein et al. (1998), Eberlein and Prause (1999)); hyperbolic distributions 
were constructed by Barndorff-Nielsen (1977)); 

• Normal Inverse Gaussian Processes (NIG) were introduced by Barndorff-Nielsen 
(1998) and used to model German stocks by Barndorff-Nielsen and Jiang (1998); 

• Truncated Levy Processes (TLP) constructed by Koponen (1995) were used for 
modeling in real financial markets by Bouchaud and Potters (1997), Cont et 
al (1997) and Matacz (2001); the extended Koponen family was constructed in 
Boyarchenko and Levendorskii (2000) (the generalization was needed since prob- 
ability distribution of Koponen's family have tails of the same rate of exponential 
decay whereas in real financial markets, the left tail is usually much fatter; in 
Carr et al (2002) and Boyarchenko and Levendorskii (2002a,b), the extended 
Koponen family is called CGMY-model and KoBoL family, respectively). 

• Normal Tempered Stable Levy processes were constructed in Barndorff-Nielsen 
and Levendorskii (2001) and Barndorff-Nielsen and Shephard (2001a); they con- 
tain NIG subclass. 

In Boyarchenko and Levendorskii (2000), a general class of Levy processes, which con- 
tained all the classes listed above modulo certain reservation about VGP was introduced, 
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under the name Generalized Truncated Levy Processes. Later, in Barndorff-Nielsen and 
Levendorskii (2001), the name: "Regular Levy processes of exponential type" (RLPE) 
was suggested. For a more detailed exposition, see Boyarchenko and Levendorskii (2002a, 
2002b). In order to present examples, recall that a Levy process can be completely spec- 
ified by its characteristic exponent, ip, definable from the equality E[e L ^ ,x<yt ^} = e~*^®. 
The characteristic exponent is given by the Levy-Khintchine formula 
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(2.1) V(0 = + 7;(M,t)+ / (l + i(^y)li.i<i(y)-e^)F(dy), 

where A := SS T is the variance-covariance matrix of the Gaussian component, b G R ri 
and F(dx) is the Levy density (density of jumps), which satisfies 



/ min{|x| 2 , l}F(dx) < oo. 

J~R n 



Any generating triplet A, b, F(dx) with these properties defines a Levy process (see e.g. 
Sato (1999)). If S = 0, then we have a pure jump process. 

Wide families of jump-diffusion processes are subclasses of the class of RLPE. In the 
first example, we introduce the family which is widely used in affme term structure 
models (see Duffie et al (2000) and Chacko and Das (2002)). 

Example 2.1. Let X be a Levy process with the Levy density 

F(dx) = c + A + e A+x T(_ OOj0 )(x)rfx + c_(— \-)e x ~ x l( 0t+oo )(x)dx, 
where A + > 0, A_ < — 1 and c± > 0. Then 

V>(£) = — £ - ib£ + - — — + 



2 s s A+ + ^ X- + iC 

where a 2 > and b G R are the variance and drift of the Gaussian component. The 
if)(£) is analytic in the strip G (A_, A+). 

Example 2.2. The characteristic exponent of a process of KoBoL family in ID is of 
the form 

(2.2) ^(0 = -itf + cT(-iz) [A+ - (A + + i£Y + (-A_) y - (-A_ - ^)1, 

where v G (0,2), u ^ l,c > 0, A_ < < A + , and /x e R; it is analytic in a strip 
3^ G (A_, A + ), and (|3.8|l - (|3.9jl are satisfied in this strip. 

Example 2.3. The characteristic exponent of a Normal Inverse Gaussian process in ID 
is of the form 

(2.3) V(0 = + ^[(« 2 - (P + ^) 2 ) 1/2 - (« 2 - P 2 ) lf % 

where v G (0, 2), 5 > 0, and a > |/3|; it is analytic in the strip G (—a +/3, a + /3), and 
f)3.8|l - (|3.9jl are satisfied in this strip, with z/ = 1. 

Since the sum of the characteristic exponents of two RLPE's is the characteristic 
exponent of an RLPE, the list of model examples can easily be expanded. For multi- 
dimensional examples, see Boyarchenko and Levendorskii (2002b). 
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Examples 2.1-2.3 are examples of pseudo-diffusions if A+, |A_|, and a ± (5 are large. 
Typically, processes observed in empirical studies of financial markets (hyperbolic pro- 
cesses and variance gamma processes including) enjoy this property. 

The majority of papers on Levy models deal with asset pricing. Eberlein and Raible 
(1999) consider the HJM-model driven by a Levy process (see also Eberlein and Ozkan 
(2001)). For the usage of jump- diffusion processes and more general Levy processes in 
affine term structure models of interest rates, see Duffle et al. (2000, 2002), Chacko 
and Das (2002) and the bibliography therein. Barndorff-Nielsen and Shephard (2001b) 
suggested to use Levy-driven Ornstein-Uhlenbeck processes for interest rate modelling 
purposes. For the subsequent developments, see Barndorff-Nielsen et al (2002). 

3. The model 

3.1. Levy-driven QTSM. In the Gaussian QTSM, the instantaneous interest rate is 
represented as a quadratic function of the state variables, and the latter are specified as 
diffusions. We assume that under an EMM chosen by the market, the SDE of the state 
variables can be written as 

(3.1) dX(t) = 0(t) - KX(t))dt + dZ(t), 

where {Z(t)} is an n-dimensional Levy process, 9 : R n — > R is a continuous vector- 
function, and k is a constant n x n matrix, whose eigenvalues Xj satisfy the condition 

(3.2) UXj > 0. 
The interest rate is modelled as 

(3.3) r(X(t)) = i? + 2(R 1 ,X(t)) + (TX(t),X(t)), 

where Rq G R, Ri G R™ are constant scalar and vector, (•, •) is the standard inner 
product in R n , and T is a positively definite symmetric matrix. The last condition 
ensures that 

r(X(t)) = (T(X(t) + F^Rj^it) + F^Rt) + Rq - WT^R^ 2 

is semi-bounded from below. By choosing Rq,R± and T appropriately, one can ensure 
any lower bound on r(X(t)). Notice that if one wishes to price a derivative of a stock 
whose dynamics is characterized by X, then one may allow r to depend only on some 
of the factors Xj(t), say, r = r(Xx(t), . . . ,X m (t)), where m < n; in this case, in (|3.3|) . 
Ri G R m , and r is an m x m matrix. 

If Z has no jump component then the bond pricing problem reduces to a system of 
ODE (Riccati equations), which can easily be solved numerically, and in the one-factor 
case, even analytically. (In the multi-factor system of Riccati equations can be 

reduced to a linear system; for the explicit realization in the framework of the Gaussian 
QTSM, see Kim (2003)). It seems unlikely that a reasonably simple exact solution exists 
for a general Levy process but we manage to obtain an asymptotic solution if X is a 
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pseudo-diffusion, that is, the Levy density of Z decays exponentially, and the rate of 
decay is large. The leading term of the asymptotics is the price in the Gaussian model 
with the same instantaneous moments of order one and two, and the correction terms 
are polynomials in the factors with coefficients depending on the time to expiry. After 
the leading term is found, they can be calculated recursively, by using only integration 
procedures in ID. Thus, the suggested method is relatively simple (though in multi- 
factor models, the number of additional integration procedures may be rather large; it 
is important that all the integrations remain one-dimensional, even in a multi-factor 
model). In the one-factor case, the first correction term is proportional to skewness, 
and the second one depends on the skewness and kurtosis; to be more precise, the first 
correction is proportional to skewness, and the second one is the sum of two terms, one 
of which is proportional to the square of the skewness, and the other to the kurtosis. In 
many cases, the contribution of the kurtosis is small relative to the other terms; if we 
omit the last term, then the pricing formula becomes a sum of the leading term which 
looks as the price in the Gaussian model, and the correction term, which is a quadratic 
polynomial w.r.t. to skewness. 

Similar formula for the forward rate and numerical examples show that the first cor- 
rection term has a pronounced upward hump, if the skewness is negative; in the result, 
the corrected forward rate curve can be hump shaped even when the Gaussian forward 
rate curve is not, and all parameters of the model are time- independent. By changing 
the parameters, various shapes of the forward rate curve can be obtained. 

Empirical studies show that both skewness and kurtosis can be fairly large, and hence, 
the corrections to the Gaussian price quite sizable. Consider, for instance, the statistics 
for the daily change interest rates idr) from Table 1 in Das (2002). (The table presents 
descriptive statistics for the Fed Funds rate over the period January 1988 to December 
1997, and the unit is 1 percent). Mean: m = —0.0005; standard deviation: a = 0.2899; 
skewness: A 3 = 0.3950; excess kurtosis: k 4 = 19.8667. Recall that for probability 
distribution P(dx), 

/+oo 
xP(dx), a 2 := ((x — m) 2 ), 
-oo 

A : = ((x - m) 3 )/a 3 , k A := ((x - m) 4 )/a 4 - 3, 

and that if P(dx) = P&t(dx) is the probability distribution of a Levy process with the 
characteristic exponent tp, then 

m(At)/At = #'(0); a 2 {At)/At = ^"(0); 

((x - m) 3 )(At)/At = -^ (3) (0); [{(x - m) 4 )(At) - 3a 4 (At)]/At = -^(0). 

We see that the coefficients in the third and fourth terms in the Taylor series for if) 
around zero are smaller than the second one but non-negligible whereas in the Gaussian 
case all coefficients starting from the third one are zero. 
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The skewness and kurtosis of the process under an EMM can assume essentially ar- 
bitrary values provided they are small w.r.t. variance; in particular, one should expect 
that the skewness of the process under EMM is negative even when the one under his- 
toric measure is positive as in the empirical example above. This means that even the 
one-factor approximate non-Gaussian model has two free additional parameters (albeit 
small) which can be used to get a better fit to the data than in the Gaussian model. In 
multi-factor models, the number of additional free parameters is larger still. 

3.2. Reduction of a pricing problem to a boundary problem. Consider a con- 
tingent claim with the maturity date T and payoff g(X(T)). Its price at time t < T is 
given by 

(3.4) f(X(t),t) = E t [ / 



exp(W r(X(s))ds) g(X(T)) 



(We consider the pricing under a risk-neutral measure chosen by the market). In appli- 
cations, the payoff g is measurable (usually, continuous), and it may grow at infinity. In 
the latter case, additional conditions on Z may be needed. For instance, if g grows not 
faster than an exponential: 

(3.5) \g(x)\<Ce"M, 

where C and u > are independent of x, then it suffices to assume that there exists 
A > uj such that for all p in the ball < A, and some t > 0, 

(3.6) E[e^ z ^\ < oo, 

which implies that the tails of probability densities of the process Z decay exponentially: 
faster than exp(— p\x\), for any p < A. 

It follows from (|3.6|) . that for any £ = 77 + it G C n from the tube domain R n + iU\ := 
{£ I = |t| < A} (in the one-factor case, a tube domain is a strip), and any t > 0, 

(3.7) E[e l ^ z{t)) ] < 00. 

(Instead of balls U\, one can use more general open sets containing the origin.) It 
is immediate from ()3.7|) . that ^(O an d its derivatives w.r.t. the complex argument £ 
are well-defined in the same tube domain R n + iU\ (one says that is analytic in 
R n + iU\), and we may use the latter condition on ip instead of the former condition 
(J3.7j) . To justify the use of the Feynman-Kac formula, we assume that Z is a regular 
Levy process of exponential type (RLPE). This means that ip admits a representation 

(3.8) m = -i(», 0+0(0, 

where p € R", and <fi satisfies the following condition: there exist c > 0, v e (0, 2] and 
V\ < v such that as £ —>■ 00 in the tube domain R" + iU\, 

(3.9) #o = cier+o(ier) 

(see Boyarchenko and Levendorskii (2002b)). The v and U\ are called the order and 
type of the process. 
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To simplify the justification of the use of the Feynman-Kac formula, we add unnec- 
essary condition: for any multi-index a = (cti, . . . , a n ), there exists a constant C a such 
that for all £ in the tube domain R ra + iU\, 

(3.10) \d a <p(o\ <c a (i + Kir' a i 

where \a\ = a± + ■ • ■ + a n . Notice that this condition holds for all model classes of 
RPPE's. 

In the appendix, by making use of the Feynman-Kac formula, we will prove the fol- 
lowing theorem. 

The orem 3.1. Let the stochastic factor satisfy / TO) . £P) . / TO) . / TO) . / TO) . 

and let r be given by \3. 3]) . and let g be a continuous function, which admits a 

bound ( TO)) . 

Then a) the stochastic expression Jff.^l ) defines a continuous function f , which admits 
an estimate 

(3.11) \f(x,t)\ <C ie "W, 

where C\ is independent of x and t < T; 

b) f is a unique solution to the following problem: 

(3.12) (d t + {B{t)-KX,d x )+L-r(x))f(x,t) = 0, t < T, 

(3.13) f(x,T) = g(x), 
where L is the infinitesimal generator of Z . 

Recall that the infinitesimal generator of the Levy process Z, L, can be represented 
in the form of a pseudo-differential operator (PDO) with the symbol —ip: L = —ip(D x ). 
A PDO A = a(D) with the symbol a acts on sufficiently regular functions as follows: 

(Au)(x) = (2n)- n [ e^a(0u(0d^ 
where u is the Fourier transform of u: 

u(0 = / e- iM u(x)dx. 

In particular, the partial derivative d x is the PDO with the symbol it;. 

3.3. Asymptotic pricing. The asymptotic pricing formulas will be derived under the 
following conditions. Assume that the characteristic exponent of the driving Levy process 
depends on a small parameter e > 0: = VK e >0 an d satisfies the following three 

conditions. First, we require that the A in the definition of the tube domain R ra + iU\ 
satisfies A >> e -1 / 2 . The next two conditions are formulated for £ in the tube domain 
R n + iU x : 

1) in the region |£| > e -1 / 2 , V ; ( e )0 admits an estimate 

(3.14) ^( e ,o>cier, 
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where v G (0, 2] and c > are independent of (e, £) in the region; 

2) in the region |£| < e -1 / 2 , tp(e,£) admits an asymptotic expansion: in the one-factor 

case, 

2 00 

(3.15) = "V* + y£ 2 - ^V -2 *, ' 

3=3 

where the coefficients are uniformly bounded: 

(3.16) \kj\< Ca 2 /2, 

where C is independent of j; in the multi-variate case, ()3.15|) is replaced with 

(3.17) = -ifoO + g 1 1^1 1' 

3=3 

where fcj(£) is a homogeneous polynomial of order j, which admits a bound 

(3.18) \k 3 ((j: T )-H0\<c\^, 

where C is independent of j. 

The asymptotic solution will be found in the following sections. Here we explain the 
main idea in the one-factor case. We look for the solution in the form 

(3.19) f = fo + eh + e 2 f 2 + --- . 
From (|3.15j) . we can formally write 

2 00 

(3.20) L = fid x + ^dl + Y, t j -%di, 

3=3 

and by substituting ()3.19j) and ()3.2()|) into ()3.12|) . we obtain a formal equality 

(3.21) (l + J2 e%- J + J2 e l f)j = 0, t <T, 



where 



a 2 

L = d t + (6(t) kx)O x + —d 2 x - r(x) 



is of the same form as the operator in the Gaussian model, and 

L l = k l+2 d l x +2 , 1 = 1,2,.... 

By multiplying out in ()3.21|) and gathering terms of the same order in e, we obtain the 
following series of problems. The leading term of the asymptotics is found from 

L f (x,t) = 0,t<T; 
fo(x,T) = g(x), 



PSEUDO-DIFFUSIONS AND QTSMS 



11 



which is the pricing problem in the Gaussian model; and the following terms are found 
step by step, by solving problems 

i 

L fi(x,t) = -Y,k j+ 2di +2 fi-j(x,t), t<T, 

3=1 

fi(x,T) = 0, 

for I — 1, 2, We believe that for practical purposes, it suffices to use an approximate 

formula ()3.19j) with terms up to order 2; this allows one to take into account the skewness 
and kurtosis. This approximate solution can be written as 

(3.22) / « / + ehh + (ek 3 ) 2 f 21 + e 2 k 4 f 22 , 

where fi, f 2 ± and f 22 solve equations 

L fi(x,t) = -d 3 J (x,t), 



L f22(x,t) = -d£f (x,t) 



L f 2 i(x,t) = -d A x h{x,t), 

54 - 



in the half-space t < T, subject to zero boundary condition. The explicit formulas for the 
bond price can be found in Section 0] and Section |5J Formula (|3.22j) may seem somewhat 
inconvenient for practical applications since it depends on the small parameter e, which 
is not explicitly specified. Notice, however, that 

e k 3 = -^( 3 )(0)/3!, e 2 k 4 = -^ (4) (0)/4!, 

and the derivatives of the characteristic exponent at can be inferred from empirical 
data - see Introduction. Thus, we may write (J3.22|) without e: 



(3.23) ^-^-(WV^/ 



22- 



By using (j3.23J) . the influence of the moments of order 3 and 4 on the price can be 
explicitly analyzed; and this influence is highly non-linear in (x, t), since the functions 
in (13~23J) are. 

If P := / is the bond price, then we can derive from (J3.23j) similar approximate 
formulas for the yield and forward rate. 

The final remark is: in order to find a current term of the asymptotics, we have to 
differentiate the previous terms, therefore an asymptotic solution with several terms 
may produce serious errors in the neighborhood of a point where the pay-off g is not 
sufficiently smooth. Indeed, one can hardly expect that a formula which is polynomial 
in x, can give a high order approximation in this case. Hence, in a neighborhood of such 
a point, a different asymptotic formula should be written: see Section |H1 
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4. Bond pricing: Gaussian model, one-factor case 

In this section, we recall the solution of problem ()3.12j) - (j3.13j) in the one-factor Gauss- 
ian case, when = — z/i£ + <r 2 £ 2 /2, and L = [id x + ^-<9 2 , and indicate the properties 
of the solution which are crucial for the asymptotic method to work. 

We assume that the interest rate is a quadratic function of the stochastic factor X(t): 



(4.1) r(X(t)) = R + 2R 1 X(t)+X(t) 2 . 

The dynamics of the stochastic factor, X, is governed by (jrS.lj) . where 9 is a scalar 
function, and k is a positive scalar. The bond price is given by (|3.4|) with g(x) = 1, 
hence it is a bounded solution to problem (j3.12J) - (j3.13|) with g(x) = 1 in the RHS of 
(EH- 

Set r = T — t, 6{t) = 6{T — t) + /i, and with some abuse of notation, write f(x, r) 
instead of f(x,T — r). We look for the bounded solution to the problem 

2 

(4.2) K + (%)-^ + y^-r(x))/(x,r) = 0, r > 0, 

(4.3) /(x,0) = g(x) 
in the form 

(4.4) f(x,r) =exp$ (z,T), 
where 

(4.5) $o(z,r) = A(t)x 2 + 5(r)x + C(r). 
By substituting ()4.4|) into ()4.2|) . we obtain 

(4.6) (exp(-$ )£ exp $ ) (x, r) - r (x) = 0, 
where 

2 

(4.7) C = -d T + (9(r) - n)d x + ya 2 , 

subject to $0(^5 0) = 0. Straightforward calculations yield the following system of ODE 
with zero initial data: 

(4.8) - A'{t) - 2kA(t) + 2a 2 A{r) 2 -1 = 0, 

(4.9) -B'(t) - kB{t) + 2ct 2 A(t)B(t) + 26(t)A(t) - 2R X = 0, 

(4.10) - C '(T)+a 2 A(r) + ^-B(T) 2 + 9(r)B(r)~R = 0. 
Equation ()4.8j) is solved by separation of variables: 

(4.") A ^ = A ^A^T^- 
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where A\ < < A 2 are roots of the quadratic equation 2a 2 A 2 — 2kA — 1 = 0, and 
to = 2a 2 (Ai — A 2 ) < 0. A(t) having being found, we can calculate B(r) from the linear 
equation ()4.9|) : 



(4.12) B(r) 
where u)\ = 2a 2 A\ — k, and 

h(r) 

h(r) 



2e^{A 2 I l {T)-A 1 I 2 {r)) 
(A 2 -A 1 )(A 2 -A l e<^) ■ 



(A 1 6(s) - R^e-^ds, 
(A 2 9(s) - i?i)e ( ™ )s ds. 



If 9 is independent of r, then /, can be calculated explicitly: 



h(r) 
h{r) 



A x 9 - R x 
A 2 6 - gi 

UJ\ — OJ 



(1-e-— ), 
(1 -e ( ^ l)T ), 



and therefore 
(4.13) B(r) 



(A 2 -A 1 )(A 2 -A 1 e^ 



A 2 (A 1 9 - R ± 



A 1 (A 2 6-R 1 ) 



jjir _ ft -±L( e ^ - e "T\ 



Finally, we find C(r) from ()4.1()j) by integration: 



(4.14) 



T 



a 2 A(s) 



a 



B{sf + 6(s)B(s) - Rq 



ds. 



If 9 is constant, then C(r) can be calculated explicitly. In order that B(t) and C(r) can 
be calculated explicitly, 9 need not to be a constant; for instance, one can use exponential 
polynomials. 

To end this section, we make the crucial remark on the properties of the solution. 
First, from (|4~TTj) . 



(4.15) 

and as r — > 0, 
(4.16) 



A{t) < 0, V r > 0, 
A(t) ~ A X A 2 - — — —t = A x A 2 2o 1 t = -t. 



'A 2 -A x 



Hence, for any r G (0, T], f(x,r) decays as exp(— tx 2 ), as x — > ±00, and /(£,t), the 
Fourier transform of /(x, r) w.r.t. the first argument, decays as r -1 / 2 exp(— £ 2 /(4t)), as 
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£ — > ±00. To be more specific, 



/+00 
e - lxi+A{ rW+B{r) X+ C(r) dx 
-00 

= -^yexp[C(r) + (£ + ^(r)) 2 /(4A(r))]. 

We conclude that for any N, in the region r G (0, T], |£| > e -1 / 2 , 
(4.18) l^/U^i^^e-^, 

where C^r is independent of r and £. Notice that the RHS of (|4.18|) is negligible. 

5. Bond pricing: Levy model, one-factor case 

5.1. The leading term of the asymptotics. We assume that (|3.14j) - (j3.16|) hold. Take 
\i and a 2 from (J3.15|) . and denote by /o the solution to the Gaussian bond pricing problem 
(j4~2>(jPfl : it is given by (jOjl . (l4~TTj) . flUPD and (@~HD, and satisfies (|4~TH|1 . flOA . (l4~T7l) 
and (|4.18jl . Introduce f 1 := f — fa. Since f and / are solutions to problems ()4.2j) - (j4.3|) 
and ()3.12|) - ()3.13|) . respectively, and ()3.15|) holds, we conclude that f 1 is the solution to 
the following problem: in the half-plane r > 0, 

(5.1) (-d T + (0(r) - kx)8 x - ij(D x ) - r(x))f\x, r) = -V^e, D x )f (x, r), 
where 

2 

V 1 {e,£) :=- ^O + y?-^ 

00 

00 

We also have the initial condition 

(5.2) / 1 (x,0) = 0. 

From (|3.16|) and (|4.18|) . the following estimate for the RHS in (J5.1|) follows: 

(5.3) 1^(6,0/0(^)1 < C oe r- 1 / 2 exp(-£ 2 /(8r)), 

where C is independent of e G (0,1) and r G (0, T]. By making the inverse Fourier 
transform, we obtain 

(5.4) \\T>i{e,D x )f {x,T)\\ C (Rxi ,T\) < Ce, 
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where C is independent of e G (0, 1). By applying the Feynman-Kac theorem to (|5.1j) - 
(15. 2 j) . the representation of f 1 in the form of the stochastic integral results: 



f 1 {x,r)=E. 



J° exp (- j r(X(s'))ds') V(e, D x )f (x, s)ds 



and from (|5.4|) . we derive an estimate 

(5.5) l/W)| <Ce, 
where C is independent of e G (0,1), a; 6 R and r G (0, T]. 

5.2. First correction term. Estimate ()5.5|) shows that /o is indeed the leading term 
of the asymptotics of / as e — > 0, and in view of ()3.15|) . it is natural to look for the first 
correction term in the form efi, where fi is the solution to the following problem: 

2 

(5.6) {-d T + (0(r) - nx)d x + °-dl - r{x))h(x,r) = -k 3 d 3 J (x,r), 

in the half-plane r > 0, subject to 

(5.7) /iM)=0. 
We look for f\ in the form 

(5.8) h(x,r) = k 3 f (x,r)h(x,r) = fc 3 e* o(!e ' T) /i(x,r), 

where <£>o is given by (J4.5j) . By substituting into (|5.6j) . we obtain that /i solves the 
problem 



(5.9) (£ + a 2 (^$ )9, + e-*°(/:e* ) - r(x))f x = -e-*°%e* , 



in the half-plane r > 0, subject to 



(5.10) /i(s,r) = 0, 

where C is defined by (|4.7|) . Equation (|4.6|) allows us to simplify (|5.9|) : 

(£ + ct 2 (2A(t)x + B{r))d x )h = -e"* ^* . 
We calculate the operator 



3 



e -*o^e*o = ( e -*°d x e 



= (d x + 2A(r)x + B(r)Y, 
and by using d x l = 0, rewrite (|5.9|l as 
(5.11) £i/i(x,r) = </ (e>t), 



16 S. LEVENDORSKII 

where 

g (x,r) : = 8A 3 x 3 + 12A 2 Bx 2 + {VIA 2 + 6AB 2 )x + 6AB + B 3 , 

2 

& := d T + (-9 1 (r) + K 1 (r)x)d x ~^-d 2 x , 

9 1 (r) := 9(r) + a 2 B(r), 
Kl (r) ■= k-2ct 2 A(t), 

and A = A(t), B = B(t). Clearly, we may look for the solution to (j5.11|) in the form of 
a polynomial in x, of degree 3, with coefficients vanishing at r = 0: 

(5.12) fi(x, t) = a 13 (r)x 3 + a 12 (r)x 2 + a n (r)x + a 10 (r). 
By substituting into (j5.11J) . we obtain a system of linear ODE: 

(5.13) a' 13 + 3Kxa 13 = 8A 3 , 

(5.14) a' 12 + 2K4a 12 - 3^013 = 12A 2 B, 

(5.15) a' u + man - 26xa 12 - 3a 2 a 13 = 12A 2 + QAB 2 , 

(5.16) a' 10 - Mil - o- 2 a 12 = 6AB + B 3 , 
which can easily be integrated step by step. Namely, let 

k 2 (s) = / Ki(s')ds'; 
Jo 

JO 

e - 2K2 (r) f T e 2"2(*) (u A(s) 2 B(s) + 36 1 (s)a 13 (s))ds, 
Jo 

e~ K2{T) [ T e K2(s) (12 A(s) 2 + QA(s)B(s) 2 
Jo 

+29 1 (s)a 12 (s) + 3a 2 a 13 (s))ds, 
(QA(s)B(s) + B(s) 3 + 9 1 (s)a 11 (s) + a 2 a 12 (s))ds. 



(5.17) 




then 




(5.18) 


ai 3 (r) = 


(5.19) 


ai 2 (r) = 


(5.20) 


on(r) = 


(5.21) 


aio(r) = 



Formulas ()5.8|) . (|5.12|l and ()5.18|1 - ()5.21|1 give the first order approximation 

(5.22) /~/o-(l + eWi) 

to the bond price. The proof similar to the one of estimate ()5.5|) albeit more involved 
shows that the error of approximation ()5.22|) is 

(5.23) |/(x, r) - f (x, r)(l + ekj^x, r))| < Ce 2 (l + \x\ 2 f 2 , 
where C is independent of e G (0,1), i 6 R and r G (0, T}. 
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Unlike (|5.5|) . we have a polynomially growing factor (1 + |a;| 2 ) 3 / 2 in the RHS of ()5.23j) . 
Notice, however, that for practical purposes, one needs to know the bond price for small 
values of r(X(t)), hence for small values of X(t), and therefore the polynomially growing 
factor (1 + |x| 2 ) 3 / 2 does not matter much. 

5.3. First correction term II: the derivation based on the change of variables. 

To simplify the calculation of the next terms of the asymptotics, it is advantageous to 
change the variables in equations similar to (|5.11|) : 

(5.24) x = -0 2 (r) + e K2(T) y, 

where k 2 is given by (J5.17j) . and 6 2 is the solution to the Cauchy problem 

6' 2 (r) - k 2 (t)6 2 (t) = 9 l (r), 
02(0) = 0, 

that is, 

e 2 ( T ) = e K2(r) f e" K2(s) 0i(s)c/s. 
Jo 

The same change of the variables simplifies the calculation of fx. Introduce an operator 
S by S(f)(y,r) = f(x(y),T). Under the change of variables ()5.24j) . — d T + (6*i(r) — 
Ki{r))d x i — > — d T and d x i— > e~ K2( - T ^d y , therefore 

2 

C 2 := S-'dS = d T - %-e- 2K ^dl 

and we can rewrite (j5.11|) in the form 

(5.25) C 2 F 1 (y,T) = G (y,r), 

where Fi = Sfi, G = Sg . Clearly, G° is a polynomial in y of the same order as g : 

G Q (y, t) = G , 3 (r)y 3 + G , 2 {T)y 2 + G ,i(r)y + G , (t), 

and the coefficients Gqj can easily be calculated by using formulas for the coefficients of 
(?o or, better, independently. Under the change of variables ()5.24|) . d x + 2A{r)x + B{j) 
becomes 

V:=e-^ T \d y + A 1 {T)y + B 1 {r)), 

where 

A x {r) = 2e 2 ^A(r), B x (t) = e^ T \B(r) - 2A{t)0 2 {t)), 

therefore 

(5.26) G = V 3 -l 

= e-^(A\y z + 3AlB lV 2 + (3 A 2 + ZA^Dy + 3A.B, + Bf), 

where k 2 = k 2 (t),Ai = A\{t) and B\ = B\{t). The solution of ()5.25|) subject to 
Fx (y, 0) = is a polynomial in y of the same order as G : 

Fi(y, t) = F h3 (r)y 3 + F h2 (r)y 2 + F hl (r)y + F 1>0 (r), 
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whose coefficients can easily be found by integration: 

(5.27) F 1)3 (r) = / Go, 3 (s)ds, 

Jo 

(5.28) F 1)2 (r) = / G , 2 (s)ds, 

Jo 

(5.29) Fx^r) = [ {G^{s) + 2,a 2 e 2 ^G^{s))ds, 

Jo 

(5.30) F 1)0 (r) = [ T (G 0fi (s) + <j 2 e 2K ^G 0i2 (s))ds. 

Jo 

After that we make the inverse change of variables y = e~ K ' 2 ^ T \x + 02 (t)), and calculate 

A = s-'f,. 

5.4. Next terms of the asymptotics. Suppose that the approximation of order j > 1 
has been found: 

/ = /<>• i+E^s/i ^/ot^/i, 



2=1 / ;=o 



where A; 2 = 1, /o = e* , / = 1, and fi, 1 < I < j, are polynomials in x with coefficients 
depending on on r: 



mi 

s 



(5.31) f l (x,r) = J2^is(r)x 



s=0 



We look for the next term of the asymptotics in the form ^ +1 kj +3 f fj + i, where : = 
fofj+i is the solution to the problem 

(5.32) (jC-r(x))f j+1 (x : r) = - 9j (x, r), r > 0, 

(5.33) f j+ i(x,0) = 0, 

where ^ +1 gj is the collection of terms of order e J+1 in the expression 

00 ( J ~\ 

p=3 \ Z=0 / 

that is, 

£?' = ^ k p k j+5 -pdP(fof j+ 3- p ). 

p=3 
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Set 

(5.34) ~g 3 = e-^ 9j 

i+3 

p=3 

i+3 

(d x + 2A(r)x + B(r)) p f 3+ ^ p , 

p=3 

and multiply (|5.32|) by e - * . We obtain 

(5.35) (£ + a 2 (^$o)^ + e~^(Ce^) - r(x))f j+1 = -~g r 

Equation ()4.6|) allows us to simplify ()5.35|) and obtain a problem in the half-space r > 
with the unknown fj+\. 

(5.36) f-d T + {9 1 (T)-K 1 (T)x)d x + Y^jfj+i = -9j,r>0, 

(5.37) f j+ i(x,0) = 0. 
From ()5.34|) . g 3 is a polynomial in x as well: 

(5.38) ~g 3 (x,r) = J2hs(r)x s , 

s=0 

where m'- = maxo<i<j(mi —l+j + 3), and any of the coefficients b 3>s may be zero, that 
is, the order of g 3 may be less than m'j. Denote by m 3 the order of g 3 . 
By making the change of variables (|5.24j) . we simplify (|5.36j) : 

(5.39) C 2 F 3+l = G 3 , 

where F 3 = Sfj, G 3 = Sg 3 . Since G 3 is a polynomial in y of order m 3 : 

Til j 

G 3 {y,T) = ^G 3 , l (T)y l 

1=0 

(the coefficients Gjj will be calculated in Subsection 15. 5j) . F 3+ i also is: 

rrij 

F j+1 (y,r) = J2 F i+^y l > 
1=0 

and the coefficients are easily found by integration: 

(5.40) F 3+1i i(t) = / G 3 j(s)ds, I = m 3 ,m 3 - 1; 

(5.41) F, +U (r) = jT (ye 2,(a «(Z + m + 1)Fj+iM*) + G ^?j da, 
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for I = rrij — 2, rrij — 3, . . . 0. 

5.5. Calculation of Gjj. We can rewrite ()5.34|) as 

3+3 

(5.42) Gj = kpk j+5 - p V p F j+3 _ p , j = 0,1,.... 

p=3 

By using the initial data F — 1, we find Gj and Fj + x from ()5.42|) and (j5.40|) - (j5.41|) step 
by step. In particular, G is given by ()5.26|) . and 

(5.43) d = k 2 3 G n + hG 12 , 
where 

(5.44) Gii = V 3 F U and G 12 = V 4 F = £> 4 • 1 
are polynomials of degree 6 and 4, respectively. Hence, 

(5.45) F x = k 2 3 F 21 + hF 22 , 

where F 2 \ and F 22 are polynomials of degree 6 and 4, respectively, which solve ()5.39j) 
with Gn and G\ 2 in the RHS (and satisfy the initial condition F 2 j(y,0) = 0). The 
reader can use ()5.40|) - ()5.41|) to obtain explicit formulas for coefficients of F 2 j in terms of 
coefficients of Gy. Notice that though the latter can be written explicitly, in practical 
implementation of the method, it is simpler to write a program which calculates the co- 
efficients of a polynomial VP, given coefficients of a polynomial P, and use this program 
to calculate G\j (and Gi for / > 1 should one wish it, though in applications, it seems 
not a reasonable thing to do). 

5.6. Second order approximation for the bond price, yield and forward rate. 

We see that the second order approximation can be written in the form 

(5.46) / w e*°(l + ekjt + (ek 3 ) 2 f 21 + e 2 kj 22 ), 

where $o, fi = S~ 1 Fi, / 2 i = 5'~ 1 -^2i and f 22 = S^ 1 F 22 are polynomials in x with 
coefficients depending on r - and parameters 8, k, fi := iip'(0), a 2 := ^"(0). For practical 
applications. ()5.46j) can be rewritten in the form 

(5.47) / « e*° (l + K 3 f\ + K 2 J 21 + Kj 22 ) , 
where 

K 3 : = eA; 3 = -# (3) (0)/3!, 
K 4 : = e 2 h = -^ (4) (0)/4! 

can be inferred from the data. 

Denote by P := / is the price of the bond, and by P := ex P(^o) the price in the 
Gaussian model; then ()5.47|) becomes 

(5.48) P w P + K 3 P X + K 2 P 21 + K 4 P 22 , 
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where 

Pi — Pofu P21 — P0J21, P22 — -£0/22- 
By using the formulas for the yield 



R(x,t) 

and forward rate 



In P(x, t) 



F(x,t) = ~\nP{x,r) 



dr 



(5.49) R^R + K 3 Ri + K\R 21 + K 4 R 



we obtain approximate formulas 
(5.49) 
where 



22- 



R (x,t) = -$ (i,r)/r, 

R x (x,r) = -fi(x,r)/r, 

R*i{x,t) = (/i(x,r) 2 /2-/ 21 (x,r))/r, 

#22^, t) = ~fzi(x,T)/T, 



and 



(5.50) F sa F + + K'F 2l + K 4 F 



22- 



where 



d 

F (x,t) = -— $o(x,t), 
d ~ 

Fi(x,t) = -—fi(x,T), 
F 21 (x,t) = A(/ 1 ( X|T )a/ 2 -/ 21 (x,r)) > 
^22(2:, r) = -—f 2 2(x,T). 

5.7. Numerical examples. (The author thanks Nina Boyarchenko for writing the pro- 
grams for numerical examples and checking the algebra in the previous two sections.) 
We take the simplest model for r: r = x 2 , and constant 9(t) = 0.06, k = 0.3, fi = and 
a 2 = 0.08. We also fix x = 0.25, and study shapes of correction terms to the bond price, 
yield and forward rates in (|5.48|) . (|5.49|) and (|5.50|) (see Figures 1-3). 

From Figure 3, we clearly see that it is the first correction term, that can account 
for the hump of the forward rate curve - provided the skewness is negative and not 
too small in modulus. In the next three Figures 4-6, we plot the bond price, yield and 
forward rate; first, the leading term (dots), then the formula with the first correction term 
taken into account (solid line), and finally, the formula with the two correction terms 
(dotted line). We take the same parameters as above, and K 3 = —a 2 / 16 = —0.005, 
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K& = K 3 /20 = 2.5- 10 -4 . We see that fairly large skewness does produce a hump-shaped 
forward rate curve, when the Gaussian curve has no hump; the asymptotic formulas are 
applicable since K 3 is small w.r.t. a 2 , and K4 is small w.r.t. K 3 . 

In the last series of figures (Figures 7-9), we fix small K4 = <t 2 /200, and show how 
the shapes of the curves vary with the skewness. 

Notice that if one fits the Gaussian QTSM and the non-Gaussian one to the same data 
set, then the leading Gaussian-like leading term in the Levy model will be determined by 
slightly different drift and volatility parameters and the spot value of the factor, hence the 
pictures above do not describe quite accurately the difference between the Gaussian and 
Levy models. However, the shapes of the curves in the Gaussian model and the leading 
Gaussian approximation in the non-Gaussian models are the same (and not differ much). 
Thus, one can use the pictures to get a feeling what a difference between the Gaussian 
and non-Gaussian QTSMs can be. 

6. Extensions and ramifications 

6.1. Multi-factor case. All the constructions in the previous two sections admit mod- 
ifications for the multi-factor case. The only two differences are 

• in the Gaussian approximation, the A(r), B(t) and C(r) are matrix functions. 
The first two can be found by solving a system of linear ODE's (for detailed 
realization in the Gaussian QTSM-model, see Kim (2003)), and then C(r) is 
found by the integration; 

• the correction terms are polynomials not in one factor but in several factors. 
This leads to systems of linear ODE's whose unknowns are coefficients at factors 
x a = x^ 1 ■ ■ ■ x% n of order a — a.\ + • — V a n < rrij, where rrij depends on the step 
of the method. 

Certainly, the systems of linear ODE to be solved become significantly larger but they 
remain linear systems; therefore, the numerical solution is fairly stable. 

6.2. Pricing under historic and risk-neutral measures. Assume that under the 
historic measure, the dynamics of X t is given by 

(6.1) dX(t) = (0(t) - K X(t))dt + dZ p (t), 

where {Z F {t)} is an n-dimensional Levy process with the characteristic exponent ip F . 
We assume that ip F admits the analytic continuation into the tube domain R™ + iU F , 
where U F is an open set containing 0. To specify the interest rate dynamics under a risk- 
neutral measure, Q, we consider first the state price deflator in the form 7r t = exp(g t ), 
where q t obeys the SDE 

(6.2) dq t = -r t dt - (A, dZ F (t)} - A t. 

The vector A e R n represent the market prices of risk of the factors. Notice that for 
processes with jumps, one cannot use an arbitrary A. For the bond to be priced, it is 
necessary that A e U F , and the condition A C U F suffices for the bond and options on 
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the yield to be priced (to be more specific, any payoff which admits a polynomial bound 
with respect to factors is admissible; if a payoff grows exponentially, then additional 
restrictions on A must be imposed). 

Set — ^ P (£ + ^A) — ip p (iA). It is easy to check that under the condition A G U F , 

ipQ is the characteristic exponent of a Levy process, call it Z®. The new measure Q is 
the Esscher transform of P popular in the literature on the pricing of options on stocks 
(see e.g. Eberlein et al. (1998) and Boyarchenko and Levendorskii (2002b)). Notice that 
ipQ is analytic in the tube domain = U p — Ad {0}, and choose A = — ip F (iA). The 
straightforward calculations show that the pricing formula 

(6.3) f(x,t) = E P 

can be written as the pricing formula under the risk-neutral measure Q: 

(6.4) f(x,t) = El 

It is possible that in some situations, the specification of risk ()6.2|) is not sufficiently 
flexible, and in the Gaussian QTSM, Ahn et al. (2002) consider a more general model 
for prices of risk. This more general specification is not applicable in Levy models with 
exponentially decaying Levy densities; however, there is an additional flexibility in model 
with jumps, which may provide any number of additional degrees of freedom: one can 
use the pricing formula (|6.4j) with any Levy process Z® provided the difference of Levy 
densities of P and Q is of finite variation (this can be shown as in the case of pricing 
of derivatives on a stock; see e.g. Carr et al. (2002)). An additional restriction should 
be taken into account if we want to use the approximate formulas of Sections 4-5: the 
Levy densities under P and Q should decay exponentially, and sufficiently fast. 

Notice that contrary to the change of measure in the Gaussian model, the change 
of measure in the non-Gaussian model may lead to the changes (albeit small) in the 
instantaneous moments of order two. 

6.3. Option pricing. Let f(x,t) be the price of an European style derivative contract 
with the pay-off g(X(Ti)), where T\ < T. A typical example is a call option on the yield 
with the pay-off g(X(T\)) = max{i?(X(Ti) — K, 0}. For a fixed T\ < T, the formulas in 
the preceding section allows one to find an approximation to g(X(T\)) as a function of 
the spot values of factors x = X(t) (and of r = T — Ti). In the Gaussian model, this 
is a quadratic function, and in the non-Gaussian model, it is a fourth-order polynomial 
with small coefficients of order 3 and 4. Hence, we can calculate the two roots of the 
equation R(x) = K, which are close to the two roots in the Gaussian approximation, by 
using simple perturbation technique. Denote these roots by xf(K), and represent the 
pay-off function g(x) in the form 

g (x) = g + (x)+g'(x) 1 



-*-g(X(T)) | X(t) = x 



•■I 

exp I - / r s ds ) g(X(T)) | X(t) = x 
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where 

g + (x) = l [x ^ +OD) (x)(R(x)-K), 

g-(x) = l { _ 00tX7] (x)(R(x)-K). 

By using the perturbation technique once again, we can calculate the Fourier transforms 
of functions g^, and reduce the problem of the calculation of f(x,t) to the family of 
problems considered in Sections 4-5, with the pay-off e %x( - instead of 1. The modification 
of the asymptotic calculations is straightforward albeit lengthy. In the end, we make the 
inverse Fourier transform, and obtain an asymptotic formula for the price f(x,t), t < T\. 
For a numerical realization of this formula, one needs to choose an appropriate grid in 
the frequency domain, and for each £ from the grid, solve a number of systems of ODE's. 
The inverse FFT finishes the job. Notice that this is a variant of the standard transform 
method (see e.g. Duffie et al. (2000), Chacko and Das (2002) and the bibliography 
therein) . 

When using this scheme, one should remember that due to the non-smoothness of the 
pay-off at the money, the approximate formulas will not work well near expiry, especially 
near expiry and strike. Fortunately, near the strike, a different approximation - and much 
simpler one - can be derived. Fix xq, and introduce 

(6.5) r(x ; x) = R + R'(x )(x - x ). 

Denote by f°(xo;x,t) the solution to the affine model with the short rate modelled by 
(|6.5|) . and the same dynamics of the factors under the risk-neutral measure and pay-off 
as in the initial QTSM. 

Theorem 6.1. f(x , t) = f°(x ; x , t) + o(r) + 0(e°°) as r = T x - t -> +0 and e -> 0. 

The proof of Theorem 16. II is relatively straightforward, and it relies on the exponential 
decay of the density of jumps. Since the solution in the affine model with jumps is well- 
known and fairly simple, we get an efficient approximation to the option price in the 
non-Gaussian QTSM. Should one wish it, the correction term to the formula f(xo,t) ~ 
f°(xo', xq, t) can be obtained as the price of a derivative security (in the same affine term 
structure model), which pays the stream of dividends at rate g(x). For in-the-money 
options, one can use a constant function g(xo) instead of g(x). 

Finally, for out-of-the-money options, a simpler approximation can be derived, in the 
form 

(6.6) f(x,t) ~ C(g;x)r, 

where C(g; x) depends only on the Levy density F(dx) of the process but not on the 
Gaussian part and drift. For the case of the call option on the yield considered above, 
essentially the same calculations as in Levendorskii (2003) give for x G (x~,xf) 



(6.7) 



/+oo 
g(x + y)F(dy), 
-oo 
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and the explicit analytic expression in terms of parameters of model classes of RLPE 
processes can be derived (the modification of the calculations in Levendorskii (2003) 
made for g(x) = (e x — 1) + and g(x) = (1 — e x ) + is straightforward). 

6.4. Parameters' fitting. Far from expiry, the leading term of the yields depends on 
the instantaneous moments of order one and two, and for out-of-the-money options on 
yields, near expiry, the leading term depends only on the jump part of the process. This 
allows us to suggest the following scheme of the fitting the model to the data, under a 
risk- neutral measure. 

1. Infer parameters of the Gaussian model (including the spot values of the factors) 
from the data on yields. Here one can use the efficient method of moments as in 
Ahn et al. (2002) or an estimation method integrated with the extraction of the state 
variables (the extended-Kalman-filter-based quasi-maximum-likelihood estimate) as 
in Kim (2003), say. We regard these parameters' values as a zero-order approximation 
to the spot values of factors, drift, mean reverting and variance-covariance parameters 
of the non-Gaussian model. 

2. Given the spot values of the factors, one can infer the conditional characteristic func- 
tion of the process from empirical data as in Singleton (1999) and calculate the mo- 
ments up to order 4, or infer these moments as in Collin-Dufresne et al. (2003). 
However, in order to identify the contribution of jumps more accurately, the following 
steps seems to be reasonable. 

3. Choose a parametrized model for the Levy density (or one of the standard Levy 
models described in Section 2), and by using the prices on interest rate derivatives 
near expiry and approximate formulas near expiry, fit the parameters of the Levy 
density, and calculate the instantaneous moments of order 3 and 4. 

4. Calculate the correction terms to the yields in the Gaussian approximation by us- 
ing the asymptotic formulas of Section 5 and the zero-order approximations for the 
parameters in the Gaussian approximation, and moments of order 3 and 4. 

5. Subtract the correction terms from the data, and use the new data set in the Gaussian 
procedure to infer the corrected values of the spot factors and the first two moments. 

6. Step 3 can be repeated in order to get a corrected specification of the Levy density 
and moments of order 3 and 4. 

Notice that in order that this procedure be consistent, the resulting moments of order 3 
and 4 must be small relative to the smallest eigenvalue of the matrix of the instantaneous 
second moments. 

7. Conclusion 

We constructed a class of QTSM models with a regular Levy process of exponential 
type in place of the Gaussian one in standard QTSM. By using the Feynman-Kac formula, 
we have reduced the pricing problem for an interest rate derivative to a boundary problem 
for a pseudo-differential operator. In the case of the bond, we have found an approximate 
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solution to the boundary problem assuming that the tails of probability densities of a 
process decay sufficiently fast. The leading term of the approximate solution looks as 
in the Gaussian model (even when the underlying process has no Gaussian component), 
and the correction terms depend on skewness and kurtosis. 

Numerical examples are produced to show that by changing skewness and kurtosis, 
various shapes of the forward rate curve can be obtained. In particular, negative skewness 
can produce a hump-shaped forward rate curve even when the Gaussian curve has no 
hump: the very non-Gaussianity of the process is (one of) causes of the hump of the 
forward rate curve. Bond prices and the yield curve also change but the types of the 
shape of the curves remain essentially the same. 

We discussed possible choices of a risk-neutral measure, and indicated additional flex- 
ibility which the usage of jumps provides. A brief outline of asymptotic pricing in the 
multi-factor case, and the pricing of interest rate derivatives is given. We derived sim- 
ple asymptotic formulas for option prices near expiry, and suggested a procedure of the 
fitting of the model to the data. 

Notice that the use of pseudo-diffusions for option pricing instead of Gaussian models 
with approximately the same number of parameters is advantageous not only near expiry 
where the latter are expected to produce serious errors. If the number of parameters is 
approximately the same, then the number of factors in the former is smaller than in the 
latter, therefore the use of the inverse FFT simplifies significantly in the non-Gaussian 
model, and the numerical procedure becomes much more stable. 



A.l. Proof of Theorem 13. 1L By using the decomposition g = g + — g_, where g+(x) : = 
max.{g(x), 0} and g>_ = g + — g are non-negative, we see that it suffices to prove Theorem 
13.11 for continuous non-negative g. Fix \ £ C°°(R) such that < x( x ) — 1 f° r a ^ x -> 
x(x) — l,x < 1, x( x ) = 0,x > 2, and for any m > 0, set g m (x) := x(\ x \/ m )9( x )- 
Then g m is a continuous function with the compact support. Define f m by ()3.4j) with 
g m in the RHS. For any x, g m (x) j g(x) asm-> oo, and by the Monotone Convergence 
Theorem, f m (x) j f{x). Notice that f m — > f in the sense of generalized functions: for 
any non-negative u G C£°(R n x (0,T)), 



Below we will show that 

(i) let g be continuous and satisfy (|3.5j) : then problem ()3.12)) - (|3.13|) has a unique con- 
tinuous solution f(g;x,t), which satisfies f|3.11|) : 

(ii) in the half-space t < T, f is of the class C 2 ' 1 w.r.t. (x,t); 

(iii) f(g m ; •, •) — > f(g; •, ■) as m -> oo, in the sense of generalized functions; 

(iv) if g is a continuous function of the compact support, then f(x,t) := f(g;x,t) is 
given by (|3.4|) . 



Appendix A. Proofs of technical results 
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By (iv), f(g m ; •,•) = / m (-, •)> by (iii) and (ii), the limit in (iii) is a continuous function, 
and since we already know that f m (x,t) — > f(x,t) pointwise, we conclude that /(-, ■) = 
/(#; •, •) solves the problem ()3.12|) - ([3.13p . and finish the proof of Theorem 13.11 

It remains to prove (i)-(iv). We start with the proof of (i). Assume that k is diag- 
onalizable: there exists a matrix C such that Kc '■= C~ x kC is a diagonal matrix with 
the diagonal entries Kj (if k is not diagonalizable, an additional step is to be made - see 
the end of the proof of (i)). By making the change of variables x = Cy, we reduce to 
the case of the diagonal matrix k(— Kc)] the becomes ^c(C) := ^((C) X 0- To 
simplify the notation below (and without loss of generality), we assume that k itself is 
diagonal. In (j3.12|) - (|3.13|) . change the variables: 

t — T — t, Xj = -9 2 j{r) + e KjT yj, j = 1, . . . , n, 

where 

§ 2j (r) = e TK > [ e'^d^T - s)ds 
Jo 

is the solution to ODE 

6' 2j (t)-k6 2j (t) = 6 j (T-t), 

subject to 6*2^(0) = 0, and set 

v(y,r) = f(x,t), 

e TK = diag(e KjT ), 
ri (y,r) = r (-9 2 (r) + e TK y). 

We obtain 

(A.l) (d T + ij(e- KT D y ) + ri (y,r))v(y,r) = 0, r > 0, 

(A.2) v(y,0) = g(y). 

Notice that 6\ e C 2 ([0,T]) since 9 G C([0,T]), and n satisfies estimate 

(A.3) |^n(y,r)| < C a , s (l + \y\) 2 - la \ |«|,s = 0,l,2, 

(A.4) c \y\ 2 -C < n(y,r), 

where Co > and Co, C ayS depend on T but not on x 6 R" and r 6 [0, T}. 

Estimates (l3~%l) . (l377|) . (l3~Tn|l . and (!A~3|l - (lA~4l) allows one to apply the standard tech- 
nique of construction of the inverse to the operator of a boundary problem for PDO to 
problem (jA.lJ) - (jA.2|) . This technique is based on the construction of an appropriate parti- 
tion of unity, localization and patching of an approximate inverse from local inverses; for 
the realization of this general scheme for many classes of PDO see Levendorskii (1993). 
In op. cit., boundary value problems in L p -based spaces were considered whereas here 
we need corresponding results for C s -based spaces. This modification is straightforward: 
see e.g. the modification in Barndorff-Nielsen and Levendorskii (2001), for a different 
class of PDO. 
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In the result, we obtain that v, the continuous solution to problem fjA.l(l - (|A.2jl . which 
admits estimate (|3.11|) . exists and it is unique. Moreover, it is of the class C 2 ' 1 in the 
half-plane r > 0, and satisfies estimate 

(A.5) sup \e-^v(y,T)\ < Csup |e~" ls/| #(?/)|, 

R n x [0,T] R n 

where C depends on T, k, 9 and the constants in estimates for ip and r. By making the 
inverse changes of variables and unknowns, we obtain (i) and (ii). 

If k is not diagonalizable, then prior to the change of variables x = Cy, an additional 
change of variables Xj i— > e Pj Xj,j = 1, . . . ,n, where pj > 0, is needed. The k will be 
replaced by k — diag(pj), which generically has pairwise distinct eigenvalues and hence, 
diagonalizable; 9 will change as well, and ip{D x ) becomes ip(e~ PlT D xl , . . . , e~ pnT D Xn ). 
After that we make the same changes of variables (using the new k and 9). 

To prove (iii), we take lji G (u, A), and apply the argument above starting with g — g m 
instead of g and U\ instead of u. Since 

sup |e- Ul|a!| (</(a:) - g m (x))\ ^0 asm^oo, 

estimate (jA.5|) implies that 

sup \e- ui[x \f(x,t)- f m {x,t))\ ^0 asm^oo, 

R n x [0,T] 

which proves (iii). 

It remains to prove (iv). We have seen that for a continuous g with compact support, 
/ is continuous in the half-plane t < T, and of the class C 2 ' 1 in the open half-plane. 
Moreover, f(x,r) decays faster than as x — > oo, for any u > —X (notice that a 

continuous g of the compact support satisfies ()3.5|) for any u, and in order that the proof 
of estimate (j3.11|) remain valid, we may use any uj > —A, negative ones in particular). 
Further, r is continuous and semi-bounded from below. These conditions are more than 
sufficient for the Feynman-Kac theorem to be applicable (for instance, at this stage, we 
can repeat the proof on p. 274 in Rogers and Williams (1994)), which gives (iv). Theorem 
13. II has been proved. 
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Figure 1. Components of the asymptotics of the bond price. Parameters: 
r = x 2 = 0.0625, 6 = 0.06, k = 0.3, fx = 0, a 2 = 0.08. 
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Figure 2. Components of the asymptotics of the yield. Parameters: r = 
x 2 = 0.0625, 8 = 0.06, re = 0.3, n = 0, a 2 = 0.08. 
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Figure 3. Components of the asymptotics of forward rate. Parameters: 
r = x 2 = 0.0625, 6 = 0.06, k = 0.3, fx = 0, a 2 = 0.08. 
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Figure 4. Bond price: leading term, first and second approximations. 
Parameters: r = x 2 = 0.0625, # = 0.06, k = 0.3, Li = 0,a 2 = 0.08, K 3 = 
-0.005, K 4 = 2.5 • 10" 4 . 
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Figure 5. Yield: leading term, first and second approximations. Pa- 
rameters: r = x 2 = 0.0625,0 = 0.06, k = 0.3, A* = 0,a 2 = 0.08, K 3 = 
-0.005, X 4 = 2.5 • 10" 4 . 




Figure 6. Forward rate: leading term, first and second approximations. 
Parameters: r = x 2 = 0.0625, 8 = 0.06, k = 0.3, n = 0, a 2 = 0.08, K 3 = 
-0.005, K 4 = 2.5 • 10" 4 . 
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Figure^ 7. Bond price: dependence on skewness. Parameters: r = x 2 = 
0.0625, 6 = 0.06, k = 0.3, fj, = 0, a 2 = 0.08, K 4 = 2 • lO" 4 . Dashes: K 3 = 
5 • 10 -3 ; dotted line: K 3 = 2.5 • 10" 3 ; dots: K 3 = -2.5 • 10 -3 ; crosses: 
K 3 = — 5 • 10 -3 . The solid line is the Gaussian curve. 
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Figure^ 8. Yield: dependence on skewness. Parameters: r = x 2 = 
0.0625, # = 0.06, k = 0.3, /i = 0, a 2 = 0.08,7^ = 2 • 10" 4 . Dashes: 
K 3 = 5 ■ 10~ 3 ; dotted line: K 3 = 2.5 • 10~ 3 ; dots: K 3 = -2.5 • 10~ 3 ; 
crosses: K 3 = — 5 • 10 -3 . The solid line is the Gaussian curve. 
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Figure 9. Forward rate: dependence on skewness. Parameters: r = 
x 2 = 0.0625,6? = 0.06, k = 0.3, A* = 0,a 2 = 0.08, K A = 2 • 10" 4 . Dashes: 
K 3 = 5- 10~ 3 ; dotted line: K 3 = 2.5- 10~ 3 ; dots: K 3 = -2.5- 10~ 3 ; crosses: 
K 3 = — 5 • 10~ 3 . The solid line is the Gaussian curve. 



